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Abstract  -  The  purpose  of  the  this  paper  is  to  apply  a  new  time  domain  scattering  model  [1]  (not  limited  to  small  wave  heights, 
small  slopes,  or  plane  waves)  to  developing  seas  which  may  be  swell-contaminated  and  which  assume  a  new  bimodal  ocean  wave 
directionality  function  [2],  This  attribute  of  the  directionality  of  the  ocean  waves  is  generated  through  nonlinear  wave-wave  interaction 
in  which  energy  near  the  spectral  peak  feeds  into  both,  shorter  and  longer  wavelength  components.  The  goal  in  this  paper  is  to  use 
these  new  models  to  further  refine  current  estimations  for  the  magnitude  and  Doppler  characteristics  of  the  acoustic  scattering  from  the 
sea  surface.  It  is  hoped  that  these  revised  estimates  can  be  used  by  other  researchers  trying  to  determine  the  characteristics  of  other 
scattering  mechanisms  expected  to  be  active  in  the  near-surface  region  of  the  sea  surface. 

I.  Introduction 

The  purpose  of  the  this  paper  is  to  apply  a  new  time  domain  scattering  model  [1]  (not  limited  to  small  wave  heights,  small 
slopes,  or  plane  waves)  to  developing  seas  which  may  be  swell-contaminated  and  which  assume  a  new  bimodal  ocean  wave 
directionality  function  [2].  More  specifically,  a  bimodal  distribution  of  the  first  kind  is  generated  through  nonlinear  wave-wave 
interaction  in  which  energy  near  the  spectral  peak  feeds  into  both,  shorter  and  longer  wavelength  components.  The  goal  in  this 
paper  is  to  use  these  new  models  to  further  refine  current  estimations  for  the  magnitude  and  Doppler  characteristics  of  the 
interface  backscattering  component  of  acoustic  reverberation  from  the  sea  surface.  These  revised  estimates  for  the  contribution 
that  the  interface  scattering  makes  to  received  frequency  spectrum  can  then  be  used  by  other  researchers  as  a  new  benchmarks 
from  which  refined  estimates  can  be  derived  for  other  scattering  mechanisms  expected  to  be  active  in  the  near-surface  region  of 
the  sea  surface. 

Although  the  frequency  shifting  and  broadening  of  acoustic  signals  backscattered  from  the  sea  surface  has  been  the  subject  of 
many  studies,  important  aspects  of  the  physics  driving  this  problem  remain  unclear  due  to  a  variety  of  modeling  and  experimental 
difficulties.  The  problem  is  particularly  interesting  and  complex  in  the  low  to  moderate  frequency  regime  (less  than  about  3  kHz) 
because  air  bubbles,  injected  into  the  near-surface  region  by  breaking  waves,  may  significantly  influence  the  magnitude  and 
possibly  the  spectrum  of  the  scattered  signal  through  two  different  competing  mechanisms  [3].  Also,  a  detailed  knowledge  of  the 
ocean  wave  spectrum  at  the  time  of  an  at-sea  experiment  is  usually  lacking  so  acoustic  modelers  often  adopt  canonical 
parameterizations  for  the  ocean  wave  field.  In  some  cases,  these  descriptions  result  in  misleading  acoustic  predictions  but,  more 
generally,  they  simply  do  not  support  the  detailed  analysis  needed  to  differentiate  between  different  near-surface  scattering 
mechanisms.  Finally,  most  acoustic  measurements  involve  moving  resources.  While  the  Doppler  effects  due  to  uniform 
rectilinear  motion  of  the  source  and/or  receiver  can  be  “subtracted  out”  leaving  only  the  frequency  scattering  due  to  the  moving 
air/sea  interface  and  other  scatterers  (e.  g.,  fish,  bubble  clouds)  near  the  ocean  surface,  this  is  an  idealization  that  requires  further 
scrutiny.  For  example,  the  undulating  motion  experienced  by  a  source  or  receiver  array  as  it  is  pulled  through  the  water  will 
cause  frequency  modulation  effects  that  tend  to  smear  out  the  Doppler  spectrum.  Although  typically  not  estimated  at  all,  this 
system  induced  spreading  of  the  reverberation  frequency  spectrum  may  become  as  large  as  the  expected  Doppler  shifts  due  to 
scattering  from  the  air/sea  interface  and  obviously  cannot  simply  be  subtracted  out.  A  reasonable  approach  to  sorting  out  the 
various  components  of  the  problem  is  to  accurately  model  one  of  the  phenomena  involved  and  assign  any  discrepancy  with  the 
data  to  other  mechanisms.  When  scattering  from  fish  can  be  ruled  out,  the  surface  scattering  component  is  usually  modeled  and 
the  environmental  parameters  needed  to  drive  bubble  scattering  models  have  been  inferred.  This  seems  to  be  a  necessary 
“bootstrapping”  process  considering  the  experimental  difficulties  of  characterizing  the  dynamic  near  surface  region  at  sub¬ 
wavelength  scales. 

II.  Modeling  Mdced  Ocean  Surfaces  and  Bimodal  Directionality 

The  common  situation  encountered  in  practical  sonar  applications  is  a  “mixed  sea”  or  “swell-contaminated  sea”,  composed  of  a 
remotely  generated  wave  field  and  a  locally  generated  wav  e  system.  These  mixed  seas  are  produced  either  by  distant  storms  or 
by  rotating  winds.  A  mixed  sea  may  have  a  double  peak  in  its  frequency  spectrum.  The  study  of  wind-driven  seas,  swell,  and 
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their  interaction  has  been  the  subject  of  many  investigations  and  an  exhaustive  review  of  the  pertinent  literature  is  beyond  the 
scope  of  this  paper.  The  goal  of  this  section  is  to  set  forth  basic  criteria  for  discussing  and  modeling  deterministic  time-evolving 
realizations  that  are  representative  of  mixed  sea  states  for  the  purpose  of  simulating  the  acoustic  backscattered  response. 

A.  Mixed  Seas 

There  exists  no  universal  or  absolute  criterion  to  distinguish  sea  from  swell  and  the  terminology  itself  is  not  standardized. 
Often,  sea  is  defined  as  the  component  driven  by  the  local  wind,  and  everything  else  is  taken  as  swell.  In  analyzing  records  of 
wave  measurements,  a  separation  frequency  (  Fs )  is  usually  defined  to  facilitate  the  distinction  [4].  Wave  components  with 
frequencies  higher  than  Fs  are  assumed  to  be  largely  due  to  the  local  wind  and  wave  components  having  frequencies  lower  than 
Fs  are  associated  with  swell.  One  practical  way  to  estimate  the  separation  frequency  is  to  assume  that  the  spectral  components 
of  the  locally  wind-generated  sea  have  their  phase  velocity  less  than  or  equal  to  the  local  wind  speed.  Therefore,  the  separation 
frequency  is  the  frequency  at  which  the  phase  velocity  equals  the  local  wind  speed.  This  concept  suggests  that  the  local  sea 
cannot  contribute  to  the  lowest  frequency  (longest  wave)  part  of  the  spectrum  of  the  mixed  sea.  On  the  other  hand,  when  the 
wind  ceases  or  changes  directions,  the  high  frequency  components  of  the  spectrum  are  the  first  to  disappear.  From  numerical 
modeling,  it  appears  that  components  of  the  wave  field  having  frequencies  higher  than  about  1.  25  times  the  peak  frequency  (i.e., 
the  frequency  where  the  ocean  wave  frequency  spectrum  is  maximum)  disappear  quickly  [5]  when  the  winds  change.  Assuming 
that  the  swell  originates  from  a  fully-developed  sea  characterized  by  a  Pierson-Moskowitz  (PM)  spectrum  [6],  the  peak  frequency 
is,  Fp™k  =  0.128g/w™e//  ,  where  g  is  the  acceleration  due  to  gravity  and  u™dl  is  the  wind  speed  at  10  m  height  (of  the  previous  or 
distant  wave  system  that  will  become  the  swell).  In  this  case,  the  separation  frequency  places  a  high  frequency  limit  on 
contributions  that  the  swell  can  make  to  the  mixed  sea  wave  spectrum,  Fswell  <  1.25  •  F^k  -  g/us^el1  . 

A  review  of  the  literature  on  the  subject  of  sea  surfaces  containing  wind-driven  seas  and  swell  shows  the  following:  (a)  High 
frequencies  components  of  a  wind-driven  sea  arise  first  but  they  are  also  the  first  to  disappear  when  wind  ceases  or  turns,  (b)  The 
swell,  once  generated,  can  propagate  huge  distances  apparently  with  little  or  no  attenuation  [5].  (c)  Clear  evidence  from  field 
measurements  of  swell  interacting  with  wind-wave  development  is  lacking  [5].  Although  in  principle,  short  wind-driven  wave 
components  can  be  suppressed  by  the  presence  of  long  waves  (swell),  the  effect  is  too  small  to  be  observed  at  sea  [7].  Another 
interesting  observation  from  numerical  modeling  is  that  the  coupling  between  wind-wave  and  swell  spectra  decreases  as  the 
angular  separation  between  the  directional  spectra  increases.  For  angular  separations  larger  than  90°  there  is  no  coupling  at  any 
frequency.  In  the  case  of  a  rotating  wind,  for  rotations  less  than  90°,  the  entire  spectrum  rotates  very  rapidly  (<  1  hr).  For 
rotations  larger  than  90°  a  new  sea  develops  while  the  old  sea  gradually  decays  [5].  On  the  other  hand,  some  hybrid  wave  models 
(e.  g.  ,  GONO,  NOW  AMO)  make  the  assumption  that  the  spectral  components  of  swell  propagate  without  loss  or  interaction 
through  a  local  wind-driven  sea,  provided  they  are  more  energetic  than  the  corresponding  spectral  components  of  the  sea  [8]. 
Those  swell  components  that  are  less  energetic  than  the  corresponding  spectral  components  of  the  sea  are  completely  destroyed 
without  any  change  to  the  energy  in  sea.  Neglecting  the  sea-swell  coupling  seems  to  be  an  acceptable  approximation  provided  the 
sea’s  peak  frequency  remains  higher  than  twice  the  peak  frequency  of  the  swell. 

B.  Ocean  Wave  Directionality 

It  is  common  in  the  oceanographic  literature  to  see  the  2D  ocean  wave  spectrum  written  as  the  product  of  a  one-dimensional 
frequency,  or  in  this  paper,  wavenumber  spectrum,  S  (K) ,  and  a  directionality  function,  D  ( [K ,  , 

^{K^)  =  S{K)D{K^).  (1) 

In  recent  years,  it  has  become  more  accepted  that  the  directionality  function  of  2D  wind-driven  seas  is  generally  bimodal  and 
symmetric  about  the  wind  direction  for  frequencies  above  about  twice  the  peak  frequency.  As  summarized  by  Hwang  et  al.  [9], 
earlier  analysis  of  temporal  measurements  by  wave-gauge  arrays  of  directional  buoys  show  unimodal  directionality.  Only 
recently,  using  maximum  entropy  and  maximum  likelihood  signal  processing  methods,  have  bimodal  features  been  extracted 
from  temporal  measurements.  In  contrast,  a  standard  2D  Fourier  transform  procedure  is  sufficient  to  bring  out  the  bimodal 
features  from  3D  topographical  images.  Ewans  [10]  has  reported  observations  of  the  directional  spectrum  for  fetch-limited  seas, 
for  the  case  of  a  well-defined  200  km  fetch,  showing  a  clear  bimodal  directional  distribution.  Ewans  also  found  that  the 
nondirectional  spectrum  associated  with  these  fetch-limited  seas  conformed  closely  to  the  JONS  WAP  spectral  shape,  and  the 
fetch  dependencies  of  the  parameters  of  the  spectrum  were  quite  similar  to  those  observed  in  the  JONS  WAP  experiment.  A 
bimodal  distribution  has  been  previously  observed  from  aerial  stereo  photography  [9].  However  the  data  quality  of  these  earlier 
stereo  photographs  was  not  very  high.  As  stated  by  Wang  and  Huang  [11],  numerical  wave  models  have  shown  that  the  presence 
of  a  bimodal  directional  spreading  is  a  robust  feature  at  wavenumbers  above  the  spectral  peak.  Recently,  Huang  et  al.  [2],  [9], 
[11]  used  high  resolution  data  from  airborne  scanning  laser  ranging  system  and  GPS  (Global  Positioning  System)  technology,  to 


acquire  3D  spatial  topography  of  ocean  surface  waves  during  a  quasi-equilibrium  wind  wave  condition.  From  the  spatial  data, 
wavenumber  spectra  were  computed  directly.  They  concluded  that  for  spatial  components  of  the  wave  spectrum  less  than  about 
twice  the  peak  wavenumber,  the  distribution  is  unimodal.  At  higher  wavenumbers,  the  spectrum  becomes  bimodal,  more 
precisely,  a  bimodal  distribution  of  the  first  kind,  forming  two  lobes  symmetric  to  the  direction  of  the  dominant  wave,  which  for 
mature  seas  (long  fetch)  is  coincident  with  the  wind  direction.  On  the  other  hand,  Walsh  [12]  found  that  the  peak  of  the  spectrum 
is  not  in  the  wind  direction  until  fetch  reaches  at  least  150  km.  The  bimodal  directivity  function  given  is  expressed  as  a  function 

n 

of  horizontal  wavenumber  K  and  azimuthal  angle  0  =  a  tan  (Ky  /Kx^.  Allowing  n<(j)<-7i  and  requiring  D{K,(jf)d(j)-  \ , 

-n 

D(K, </>)=- 0.5+  an(K)cos(n<p)  .  (2) 

^  n= 1 

As  demonstrated  by  Wang  and  Hwang,  only  eight  terms  are  necessary  to  adequately  describe  the  distribution.  The  Fourier 
coefficients,  an  where  K  =  K/K peak  ,  are  calculated  using  the  approximation  an(K)  =  ClnKCln  exp  (~KC3n  C4n ,  and  the 
empirical  coefficients  Cln ,  C2n  ,  C3n  ,  and  C4n  are  provided  in  [13].  Fig.  1  shows  the  surface  roughness  spectra  for  developing 
2D  seas  that  assume  bimodal  and  the  more  familiar,  but  less  realistic,  cos2  (0/2)  azimuthal  dependence.  The  JOSWAP  spectrum 
is  assumed  as  are  winds  of  8  m/s  wind  (w10)  for  4  hrs  into  the  0  =  0°  or  Kx  =  Ky  =0  direction. 
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Fig.  1  Roughness  spectra  of  2D  developing  ocean  waves  assuming  winds  at  8  m/s  for  4  hours.  In  (a)  a  bimodal  distribution  has  been  assumed 

and  in  (b)  a  cosine-squared  azimuthal  dependence  has  been  assumed. 

A.  Generating  Deterministic  Sea  Surface  Realizations 

The  numerical  surfaces  used  in  this  study  are  generated  by  linear  filtering  in  the  2D  wavenumber  domain  an  array  of  random, 
uncorrelated  numbers  having  unit  variance  and  a  Gaussian  probability  distribution  function  (PDF)  of  surface  elevations.  The 
transfer  function  of  the  linear  filter  is  proportional  to  the  square  root  of  the  2D  power  spectrum  of  the  surface  elevations, 

Om  ( Kx ,  Ky ) ,  which  is  derived  from  the  wavenumber-angle  representation,  ( K ,  0)  using  K  =  yj Kx  +  Ky  , 

0  =  arctan  (Ky  j Kx  ) ,  and  the  Jacobian  of  the  transformation,  (1 /K).  The  time  evolution  of  the  surface  realization  is  achieved  by 
advancing  in  the  wavenumber  domain  the  phase  of  each  ( Kx,Ky )  component,  according  to  its  frequency,  known  through  the 

dispersion  relation  F  =  ±yfgK  /hr ,  and  the  desired  time  lag.  This  is  by  now  a  familiar  process  that  has  been  thoroughly 

described  in  the  literature  of  scattering  from  rough  surfaces  determined  by  stochastic  process. 

Based  on  the  considerations  identified  in  the  previous  section,  the  power  spectra  compatible  with  mixed  seas  are  constructed 
under  the  following  assumptions  and  approximations:  First,  the  wind-driven  component  of  the  surface  roughness  is  restricted  to 
developing  or  growing  seas  (duration-limited).  To  avoid  the  bimodal  directionality  of  the  second  kind,  a  phenomenon  not 
modeled  in  the  distribution  function  of  Wang  and  Hwang,  the  fetch  for  the  growing  sea  will  be  assumed  to  be  larger  than  150  km. 


The  non-directional  part  of  the  2D  spectrum  is  described  by  a  JONS  WAP  spectrum  [14],  written  here  as  a  function  of 
wavenumber: 


S,(K)  = 


a 

2 


K  3  exp 


5 

4  K 


2 

+  Di(y)exp  - 


1 

2a2 


(3) 


Here,  y=  3.3 , Km  =  (49^2^_066)/g  ,  a  =  0.076/£ 022 ,  o  =  0.07  for  K  <  Km  and  o  =  0.09  for  K>Km,  and  %=gT/us™  where 

T  is  the  duration  in  seconds  and  u^a  is  the  local  wind  speed  at  a  height  of  10  m.  Second,  the  swell  components  of  the  mixed  sea 
are  assumed  to  be  the  remnants  of  a  fully-developed  sea  generated  at  an  earlier  time  or  at  a  distant  location  and  are  represented  by 
a  truncated  Pierson-Moskowitz  spectrum  [6] : 

SPM  (K)dK  =  -  K~3  exp  -P - £ - -  dK  .  (4) 

K  1  2  F  K2  (0.92urodl) 

And  here,  a- 0.0081,  jB  =  0.74,  u™el1  is  the  wind  speed  that  generated  the  swell  measured  at  a  height  of  10  m,  and  only 
wavenumbers  Kswell  <1.01  gf  are  allowed  to  contribute.  The  angular  distribution  function  of  Hwang  and  Wang  is  also 

applied  to  the  swell  spectrum,  however,  since  the  swell  spectrum  consists  only  of  the  low  frequency  band  around  the  spectral 
peak,  this  distribution  function  reduces  to  a  unimodal  form.  Third,  the  swell  is  assumed  to  propagate  without  damping.  The 
modeling  presented  herein  cannot  be  applied  to  shallow  water  environments  where  the  ocean  waves  may  be  damped  from 
shoaling.  Fourth,  the  coupling  between  sea  and  swell  is  ignored.  Following  the  assumption  made  in  some  hybrid  wave  models, 
depending  on  their  relative  spectral  amplitudes  the  swell  either  dominates  or  disappears.  To  form  the  spectrum  of  the  mixed  sea 
Om  {Kx ,  Ky )  containing  sea  and  swell,  the  following  logic  is  applied:  At  each  point  in  the  2D  wavenumber  domain,  the  spectral 

amplitudes  <Dsefl  (Kx,Ky)  and  Q>mell  (Kx,Ky)  are  compared,  if  >®sea ,  then  d>M  =  <&„,„  ,  otherwise,  ®M  =  <&„, . 


III.  Modeling  Acoustic  Scattering  for  Moving  Ocean  Waves 

In  a  recent  paper  [1],  Keiffer  et  al.  introduced  and  benchmarked  an  approach  to  modeling  the  acoustic  Doppler  from  two- 
dimensional,  time-evolving  surfaces  that  started  with  the  time-varying  impulse  response  function.  The  method  assumes  a  point 
source  and  works  directly  in  the  time  domain  so  that  the  response  of  broadband  signals  can  be  naturally  incorporated.  It  can  also 
accommodate  realistic  source  and  receiver  motion.  The  approach  uses  time-variant  linear  filter  theory  as  a  general  framework 
and  deals  with  the  kinematics  of  the  problem  in  the  time  domain  where  the  time -varying  impulse  response  of  a  moving  interface 
is  the  fundamental  quantity  to  be  calculated.  The  impulse  response  of  the  time-evolving  surface  is  calculated  via  the  wedge 
assemblage  (WA)  model,  which  is  based  on  the  exact  wedge  diffraction  solution  of  Biot  and  Tolstoy  [15].  It  is  not  limited  by 
wave  height,  slope,  or  surface  curvature.  In  [1],  the  methodology  and  validation  for  extending  the  WA  model  to  scattering 
problems  in  which  the  surface,  source,  and  receiver  are  in  motion  was  presented.  A  systematic  study  was  performed  using  a 
series  of  increasingly  complicated  moving  surfaces.  Space  does  not  allow  for  a  complete  model  description  here.  Any  reader 
interested  in  a  time  domain  methodology  for  accurately  approximating  the  impulse  response  of  surface  or  objects  can  find  these 
details  in  Ref.  [1]. 

A.  Kinematics  and  the  Single-Scatter  Impulse  Response 

Suppose  that  a  moving  (subsonic)  point  source,  having  instantaneous  location  Rs(t\  stops  at  time  t-X  and  immediately 

begins  injecting  mass  at  a  constant  rate  into  a  homogeneous  (static)  fluid  medium  [footnote].  The  pressure  impulse  that  radiates 
outward  from  the  source  travels  (with  speed  c )  from  the  location  of  the  source  at  the  time  of  emission  and  insonifies  a  moving 
(subsonic),  impenetrable,  rough  surface,  (x,y,£).  As  the  impulsive  source  wave  front  sweeps  across  the  time  varying  surface, 

it  intercepts  or  detects  each  point  on  the  surface  at  a  discrete  instant  of  time.  This  time  of  detection,  td ,  solves  the  equation, 

td-z-\xx+yy+zs(x,y,td)z-Rs(r'\/c  =  0 .  (5) 

The  locus  of  surface  points  detected  by  the  impulsive  source  signal  defines  a  new,  effective  surface,  C(x,y,  t),  the  frozen, 
quasi-static  surface,  £ (x, y,  r)  =  (x, y,t  =  td ) . 

Causality  requires  that  the  time  history  of  the  surface  prior  to  the  arrival  of  the  incident  impulse  cannot  affect  the  scattered 
response.  Similarly,  the  single-scatter  component  of  the  scattered  response  cannot  be  affected  by  the  time  history  of  the  surface 
after  the  initial  source  insonification.  That  is,  the  single-scatter  component  of  the  scattered  response  depends  only  on  the  shape  of 


the  detected  surface.  The  effect  of  source  and  surface  motion  on  the  time  variation  of  the  impulse  response  function  can  be 
modeled  (in  the  single-scatter  approximation)  by  examining  the  impulse  responses  of  a  time  sequence  of  the  quasi-static  surfaces, 
£(x,y,ry  Receiver  motion  can  also  be  included.  Regarding  each  point  on  the  surface  £(x,y,  t)  as  a  secondary  source  that 
emits  a  scattered  impulse  at  t  =  td,  let  tR  be  the  time  at  which  that  scattered  signal  reaches  the  moving  receiver  located  at 
tR  =  td+\xx+yy+£(x9y9  t)z-Rr  (tR)\/c  .  Each  secondary  source  on  £(x,y,  r)  emits  a  scattered  impulse  that  reaches  the  moving 

receiver  at  a  discrete,  potentially  different,  future  location.  The  effect  of  receiver  motion  can  be  included  in  the  modeling  by 
using  the  static  effective  receiver  location  in  the  impulse  response  calculation  of  each  secondary  source  and  then  summing  all  the 
impulse  responses  with  due  respect  to  time  of  arrival. 

B.  Numerical  Implementation 

Although  it  is  appropriate  to  repeat  some  of  the  details  of  a  numerical  implementation  of  the  WA  model  here,  the  focus  of  this 
section  is  on  the  numerical  implementation  of  the  methodology  being  proposed  to  extend  the  WA  model  to  time-variant 
problems.  Begin  by  assuming  that  the  time  varying  surface,  ( x9y9t ) ,  is  represented  on  a  fixed  uniform  spatial  grid  (  Ax  =  Ay  ) 

and  that  the  time  evolution  of  the  surface  is  known  at  all  times.  Then,  at  any  discrete  horizontal  location,  (x/ ,  yj ) ,  the  value  of  the 

surface  height  function  detected  by  the  impulse  emitted  by  the  source  at  time  t  =  x  defines  C,  (xt ,  yj ,  r) .  Now,  the  surfaces  of 

interest  in  this  paper  are  modeled  by  filtering  random  numbers  in  the  wavenumber  domain  and  then  Fourier  transforming  to  the 
spatial  domain.  And  the  time  evolution  of  the  interface  is  accomplished  by  advancing  in  time  the  phase  angle  of  different 
(wavenumber)  spectral  components  at  rates  that  are  commensurate  with  the  gravity  wave  dispersion  relation  £l=±JgK  .  To 
numerically  determine  the  quasi-static  surface  £(xi9yj9r),  the  surface  zs(xi9yj9t^  is  first  generated  at  regular  time  intervals, 
tj  =  (J  -l)Ar  .  The  height  of  the  surface  when  the  impulse  from  the  source  arrives  is  determined  by  interpolating  between  the 
two  “snapshots”  that  bracket  the  solution  to  Eq.  5. 

To  calculate  the  (finite)  impulse  response  of  this  surface  using  the  WA  method,  the  surface  g(xi9yj9r)  is  “tiled”  using 

triangularly  shaped  facets.  These  triangular  facets  are  determined  by  first  interpolating  between  any  four  surface  heights  of  the 
square  sample  area  to  create  a  center  surface  height  value.  This  divides  the  sample  area  into  triangularly  shaped  pieces.  Adjacent 
triangular  facets  are  used  to  define  finite-length  wedge  apices.  The  impulse  response  of  the  rough  surface  is  modeled  from  the 
diffracted  responses  (P^)  of  every  (non-shadowed)  finite-length  apex  in  the  assemblage  of  wedges.  The  (finite)  impulse  response 

is  calculated  at  discrete  times,  ti  =  (i-l^At ,  using  the  exact  solution  of  Biot  and  Tolstoy  [16]  and  added  with  due  respect  to  time 
of  arrival: 

j  M  tj+At 

h  ,r)=v  foni’t’  T)h  ■  (6) 

For  a  fixed  source  emission  time  x ,  sampling  the  impulse  response  function  in  time  t  ,  at  interval  At ,  determines,  after 
Discrete  Fourier  Transform  (DFT),  a  time-variant  transfer  function,  H  (ji9  t)  ,  and  (discrete)  transform  variable,  f.  .  H  (ji9  r)  is 
defined  on  principal  interval,  -fNyq  <  ft  <  fNyq ,  the  highest  frequency  of  which  is  the  Nyquist  frequency  for  the  received 
frequency:  fNyq  =  1/2 At .  The  finite  extent  of  the  rough  surface  ultimately  limits  the  length  of  impulse  response  function  and 
determines  the  received  frequency  resolution  Af  =  l/NAt .  The  discrete  received  frequencies  are  given  by  f.  -  ~fNyq  +  (/  -  1)A/ 

where  i  is  a  positive  integer  that  is  no  greater  than  N +  1 .  Based  on  past  experience  with  static  sinusoids  and  sea  surfaces, 
choosing  At  so  that  fNyq  is  about  an  order  of  magnitude  larger  than  the  highest  received  frequency  of  interest  (/max )  is  usually 
more  than  adequate. 

Similarly,  for  a  fixed  reception  frequency,  f.  ,  sampling  the  transfer  function  H  (ji9  r)  as  it  varies  with  source  emission  time 
X  (M  samples  at  interval  Ax  )  introduces,  after  the  corresponding  DFT,  another  transform  variable  that  can  be  identified  as  the 
discrete  source  frequency,  vj=—vNyq+(J  —  V)Av .  Here  j  is  a  positive  integer  that  is  no  greater  than  M  + 1.  The  discrete 

bifrequency  function,  B  ,  Vj ) ,  that  results  after  this  second  DFT  is  specified  in  a  two-dimensional  frequency  space  over  a 
rectangular-shaped  principal  band  that  is  limited  in  one  dimension  by  -fNyq  <  ^  <  fNyq  and  the  other  by  ~VNyq  <  Vj  <  vNyq .  Here, 
the  source  Nyquist  frequency,  vNyq  =1/2 At  ,  should  be  large  enough  to  satisfactorily  minimize  the  “folding”  of  high  frequency 
information  down  into  the  principal  frequency  band. 


IV.  Numerical  Simulations 


In  this  section,  frequency  spectra  of  signals  backscattered  from  different  time-evolving,  2D  seas  are  calculated,  and  presented 
in  a  series  of  increasingly  complex  and  realistic  ocean-acoustic  simulations.  In  the  first  simulation,  the  azimuthal  dependence  of 
the  frequency  spectra  of  signals  backscattered  (30°  grazing)  from  2D  seas  that  have  assumed  the  bimodal  distribution  function  are 
examined  and,  for  reference,  compared  with  simulations  that  assumed  the  cos2  (0/2)  azimuthal  dependence.  In  other  numerical 

experiments,  the  grazing  angle  is  lowered  from  30°  to  10°,  the  effects  of  wind  duration  and  swell-contamination  are  investigated. 
In  all  of  these  numerical  experiments  (unless  specified  otherwise),  the  simulated  2D  sea  surfaces  were  specified  at  discrete  times 
on  a  finite-sized  (  Lx  =  Ly  =  251.3  m )  uniform  grid  (Ax  =  Ay  =  0.245  m ).  Each  time-evolving  surface  was  examined  for  30  s  at  0. 

2  s  intervals.  Ten  realizations  of  the  frequency  spectrum  of  the  backscattered  signal  were  computed  (one  each  from  10 
realizations  of  a  particular  sea  surface  spectrum)  and  averaged.  The  average  frequency  spectrum  is  plotted,  more  accurately, 

10 •  Logl0  (Y67T2  (rR+rs^)  (Pout  (jj)P*ut  is  plotted.  Here  7]  =  f-v  is  the  Doppler  shift,  /  is  the  received  frequency,  and  v  is 

the  source  frequency  (10  Hz  to  330  Hz).  The  factor  16  7T1  (rR  +  r5)2  comes  from  the  normalization  of  the  scattered  time  signal  by 
the  magnitude  of  the  corresponding  image  reflection. 

A.  Bimodal  versus  Cosine-Squared  Directionality 

To  examine  the  impact  that  a  bimodal  dependence  for  the  ocean  waves  may  has  on  the  acoustic  scattering,  the  azimuthal 
dependence  of  the  average  (monostatic)  backscattered  frequency  spectrum  was  calculated  for  developing  2D  seas  due  to  an  8  m/s 
wind  (m10)  that  persisted  for  4  hrs  into  the  0=0°  direction  (see  Fig.  1).  These  calculations  were  carried  out  for  the  bimodal  and 

cosine-squared  spectra  for  six  azimuthal  angles:  0°,  15°,  30°,  45°,  60°,  and  90°  (crosswind).  Note  that  an  azimuthal  angle  of 
0=0°  corresponds  to  case  of  backscattering  from  seas  whose  dominant  direction  (identical  to  the  wind  direction  for 
wavenumbers  K  <  2 Kpeak )  is  toward  the  source/receiver  location.  The  source  and  receiver  remained  fixed  in  these  experiments 
so  an  azimuthal  angle  of  0  =  90°  corresponds  to  case  of  backscattering  from  seas  whose  dominant  direction  is  perpendicular  to 

the  incident/scattered  direction.  For  azimuthal  backscatter  angles  0°,  15°,  30°,  45°,  Fig.  2  shows,  respectively,  the  average 
backscattered  frequency  spectra  (Doppler  Spectrum).  Fig  3.  shows  the  Doppler  Spectrum  for  azimuthal  orientation  of  60°,  and 
90°.  The  color  bar  in  Fig.  3  applies  to  Fig.  2  as  well.  In  all  cases,  the  nominal  grazing  angle  is  30°.  Certainly  the  most 
noticeable  differences  between  the  two  sets  of  results  are  in  the  down-shifted  parts  of  the  Doppler  spectra.  These  contributions  to 
the  received  spectrum  are  due  to  ocean  wave  components  that  are  traveling  away  from  the  source/receiver  location.  From  the 
bimodal  roughness  spectra  in  Fig.  1,  it  can  be  seen  that  for  wavenumbers  less  than  about  2.  0  rad/m,  there  are  prominent  “holes” 
in  the  Ky  <  0  part  of  the  roughness  spectrum.  It  can  also  be  noted  that  under  the  bimodal  assumption,  there  is  clear  acoustic 

evidence  in  the  0=0°  geometry  of  ocean  wave  components  traveling  against  the  wind.  In  contrast,  and  as  expected,  the  cosine- 
squared  distribution  function  yields  a  backscattered  frequency  spectrum  that  shows  no  evidence  of  ocean  waves  traveling  away 
from  the  source/receiver  location.  For  the  cos2  (0/2)  directionality  function,  the  first  evidence  of  down-shifted  contributions 
begin  with  the  0  =  15°  orientation. 

Other  features  of  the  results  worth  noting  are:  1)  The  very  weak  scattered  response  for  source  frequencies  below  about  11  Hz, 
evident  for  both  azimuthal  models  that  is  due  to  the  lack  of  ocean  wave  energy  at  the  long  wavelengths.  2)  The  spreading  about 
the  spectral  peak  is  not  symmetric.  When  there  is  a  dominant  up  or  down-shifted  Bragg-scale  ocean  wave  component,  the 
backscattered  spectrum  tends  to  spread  toward  Doppler  shifts  that  are  higher  than  the  Bragg-shift.  When  the  azimuthal 
orientation  is  such  that  there  are  comparable  Bragg-scale  ocean  wave  components  traveling  both  toward  and  away  from  the 
source/receiver  direction  (say  from  0  =75°  and  greater)  then  the  spreading  about  the  Bragg  component  tends  to  slightly  favor 
Doppler  shifts  that  are  lower  than  the  expected  Bragg-shift.  3)  For  source  frequencies  below  about  100  Hz,  it  can  be  see  that  the 
combination  of  significant  ocean  wave  energy  moving  both  toward  and  away  from  the  source/receiver  leads  to  scattering 
contributions  that  fall  into  unshifted  (  Tj  =  0  Hz)  Doppler  band. 
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Figure  2:  Azimuthal  dependence  of  monostatic  backscatter  Doppler  spectrum  from  seas  having  a  bimodal  and  cosine  -  squared  distributions. 

Grazing  angle  is  30°  and  azimuthal  angles  are:  0°,  15°,  30°,  45°. 

B.  Wind  Duration 

As  is  well  known,  when  a  (constant,  uniform)  wind  starts  blowing  on  the  ocean  surface  the  short  wavelength  components  of  the 
ocean  wave  spectrum  are  quickly  established  and  then  saturated.  As  the  wind  continues  to  blow,  longer  and  longer  wavelength 
components  become  energized,  the  peak  of  the  ocean  wave  spectrum  moves  to  lower  wavenumbers,  and  the  overall  roughness  of 
the  air  sea  interface  increases.  This  suggests  that  there  would  be  little  change  to  the  Doppler  spectrum  in  a  backscatter  geometry 
for  source  frequencies  above  several  tens  of  Hz.  Two  obvious  caveats  come  to  mind:  First,  unlike  the  cosine-squared  distribution 
function,  the  bimodal  distribution  function  depends  on  the  wavenumber  where  the  ocean  wave  spectrum  peaks.  Second,  the 
importance  of  non-Bragg  scattering  (i.  e.  ,  scattering  contributions  not  described  by  the  first  order  term  of  the  perturbation  theory) 
will  increase/decrease  with  increasing/decreasing  ocean  wave  roughness.  Simulation  results  for  the  average  backscattered 
Doppler  spectrum  (30°  grazing  and  0°  azimuth)  from  an  a  ocean  wave  field  due  to  an  8  m/s  wind  but  having  a  duration  of  only  1 
hr  are  shown  in  Figure  4a.  To  facilitate  the  comparison,  the  30°  grazing  and  0°  azimuth  calculation  belonging  to  Fig.  1  (8m/s  and 
4  hr  duration)  is  reproduced  in  Figure  4b  on  a  somewhat  larger  scale. 

Looking  at  Figure  4,  it  can  be  seen  that  when  the  duration  of  the  wind  is  shortened  the  peak  in  the  scattered  response  moves  to 
a  higher  source  frequency  (40  Hz  up  from  20  Hz  for  the  4  hour  duration  case).  The  change  in  the  bimodal  distribution  caused  by 
the  change  in  the  wind  duration  can  be  clearly  seen  by  comparing  the  down-shifted  part  of  the  Doppler  spectrum.  In  Figure  4b, 
the  absence  of  significant  long  wave  energy  makes  more  distinct  the  “turning  on”  (at  a  source  frequency  of  about  60  Hz  and  again 
at  150  Hz)  of  the  higher  order  scattering  effects  that  create  the  first  and  second  shoulders  to  the  right  (greater  Doppler  Shifts)  of 
the  up-shifted  Bragg  branch.  Note  that  these  shoulders  also  occur  in  the  plot  for  the  4  hour  duration  seas  but  emerge  at  lower 
source  frequencies  (30  and  90  Hz)  and  with  such  a  smaller  separation  in  Doppler  shift  from  the  main  Bragg  branch. 
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Figure  3:  Azimuthal  dependence  of  monostatic  backscatter  Doppler  spectrum  from  seas  having  a  bimodal  and  cosine  -  squared  distributions. 

Grazing  angle  is  30°  and  azimuthal  angles  are:  60°  and  90°  (cross-wind). 
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Figure  4:  Effect  of  wind  duration  on  Doppler  spectrum  assuming  bimodal  seas:  (a)  wind  duration  =  8  hrs,  (b)  wind  duration  =  1  hr. 


C.  Swell  Contamination 

As  stated  in  the  opening  paragraphs  of  this  paper,  the  common  situation  encountered  during  an  acoustics  experiment  is  that  of  a 
mixed  sea.  The  presence  of  swell  obviously  contributes  energy  to  the  low  wavenumber  part  of  the  ocean  wave  spectrum  that,  in 
moderate  grazing  angle  backscatter  geometries,  should  make  available  first  order  or  Bragg-scale  scattering  contributions  at 
correspondingly  very  low  source  frequencies  (perhaps  less  than  20  Hz).  The  question  more  apropos  of  the  typical  at-sea 
experiment  is  how  does  the  presence  of  the  swell  affect  the  Doppler  spectrum  at  higher  source  frequencies?  The  scenario  studied 
next  begins  with  the  wind  blowing  constantly  at  14  m/s  for  many  hours  into  the  (/)  =  180°  direction  and  then  suddenly  decreasing 
to  8  m/s  and  changes  direction  to  the  (/)  =  0°  direction.  Figure  5  shows  the  average  Doppler  spectrum  that  that  results  1  hr  after 
this  new  wind  direction  is  established.  Here  the  source/receiver  grazing  angle  is  again  30°  with  azimuthal  orientation  0=0°. 
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For  reference,  this  result  is  compared  to  Figure  4b  because  it  shares  the  same  instantaneous  wind  conditions  and  scattering 
geometry. 

From  Figures  5  it  can  be  seen  that  the  swell  contributes  to  the  acoustic  response  near  a  source  frequency  near  5  Hz.  At  slightly 
higher  source  frequencies,  between  10-20  Hz,  there  is  a  drop-out  in  the  scattered  response  because  the  swell  cuts-off  and  the 
duration  of  the  new  wind  is  too  short  to  generate  significant  ocean  waves  at  the  long  wavelengths  needed  for  Bragg  scattering  at 
10-20  Hz.  Above  a  source  frequency  of  20  Hz  the  scattering  picks  up  again  as  the  longest  ocean  waves  associated  with  the  new 
wind  direction  begin  to  supply  significant  Bragg-scale  roughness.  The  presence  of  the  long  swell  waves  increases  the  overall 
roughness  of  the  surface  and  this  causes  the  “turning  on”  of  the  higher  order  scattering  contributions  at  lower  source  frequencies. 
The  distinct  shoulders,  so  clearly  seen  in  Doppler  spectrum  shown  in  Fig.  4,  are  now  smeared  out  because  the  higher  order 
contributions  that  involve  both  the  long  waves  of  the  swell  and  Bragg  scale  components  of  the  sea  have  Doppler  shifts  that  are 
only  slightly  different  from  the  Doppler  shifts  associated  with  the  Bragg  scale  components.  Finally,  it  can  be  seen  that  the 
direction  of  the  swell  changes  the  spreading  about  the  Doppler  peak.  In  this  case,  the  spreading  is  symmetric  about  the  peak  shift. 
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Figure  5:  Effect  of  swell  contamination  on  monostatic  backscatter  Doppler  spectrum.  Swell  due  to  14  m/s  wind  added  to  seas  driven  by  8  m/s  wind  for  1  hour. 
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